function [log_likelihood,nan_pct, av_sim_prefs]=ml_risk_time_score_types_cdf(theta)
%%% calculates the log-likelihood allowing multiple factors. Joint
%%% estimation of risk and time structural coefficients



 global  X_ind_s scores ind_id factor_n  ind_n draw_n   theta_xy lottery_n lottery_choice_ind time_choice_n   time_choice_ind   n_type sim_factor_contribs maxima minima  r_scale sig_th_lower_limit CRRA2  K_time   x1 t1 x2 t2 hyper CARA sig_th_scale analytical

    
    %%% simulated factors from previous step
    sim_factors=scores;
    
    %%% Probability contribution from previous step (likelihood of measures
    %%% observed for each individual given his characteristics and
    %%% simulated draws)    
    sim_l_contrib_all_factors=sim_factor_contribs;
    
    % # of observations including rows with simulated factor draws for each
    % person
    rows=size(sim_factors,1);
    

    sex=X_ind_s(:,1);
    
 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Lottery Part %%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% All Types  %%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%% parameter assignment for risk analysis
    
    theta1=theta(1);
    theta2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];        

    kappa1=theta(1);
    kappa2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];     

    sig_th1=theta(1);
    sig_th2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];  
    
    type1_3=NaN(3,n_type);
    p_type=zeros(n_type,1);

    %%% parameter assignment for time analysis

    r1=theta(1);
    r2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];        

    sig_r1=theta(1);
    sig_r2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[]; 

    %%% For robustness with a separate mistake intercept for intertemporal
    %%% tasks
    if K_time==1
        %%% adds the extra time intercept for kappa to the end of the
        %%% coefficient list
        kappa_extra=theta(end);
        theta(end)=[];
    end    
        
    type4_5=NaN(2,n_type);

    sim_l_contrib_lottery=NaN(ind_n*draw_n,lottery_n);
    sim_l_contrib_time=NaN(ind_n,time_choice_n);
    
    av_l_part=NaN(ind_n,n_type);
    theta_i=NaN(ind_n*draw_n,n_type);
    kappa_i=NaN(ind_n*draw_n,n_type);
    
    %%% weigh over unobserved types, each has its own set of intercepts for
    %%% the structural parameters of interest
    for tt=1:n_type
        %%% last ype's proba is just 1-sum(other probas) so don't calculate
        if tt<n_type
            type1_3(:,tt)=theta(1:3);
            type4_5(:,tt)=theta(4:5);
            %%% probability needs to be between 0 and 1
            p_type(tt)=normcdf(theta(6));
            theta(1:6)=[];      
        elseif tt==n_type
            type1_3(:,tt)=theta(1:3);
            type4_5(:,tt)=theta(4:5);
            theta(1:5)=[];  
            p_type(tt)=1-sum(p_type);
        end
        
        %%% individual risk aversion based on sex and particular factor draws -
        %%% deterministic part - note that as above, the factor is determined
        %%% by the orthogonal draw AND by X*alpha - a function of observables

        %%% Even though the following parameters are individual-specific - thus
        %%% _i - they vary here within individuals due to the different factor
        %%% draws used in each row for the factor simulation
        theta_i(:,tt)=type1_3(1,tt)+sex*theta1+sim_factors*theta2_n;

        %%% individual trembling based on sex and particular factor draws -
        %%% normcdf() to bind it between 0 and 0.5 now (say that one can
        %%% make mistakes at most 50% of time, otherwise, probably
        %%% preferences misestimated and kappa takes on their role).
        kappa_i(:,tt)=0.5*normcdf(type1_3(2,tt)+sex*kappa1+sim_factors*kappa2_n);       
        kappa_i=min(maxima(2), kappa_i);
        kappa_i=max(minima(2), kappa_i);
        
        %%% individual variability of risk aversion based on sex and particular
        %%% factor draws 
        sig_th_i=sig_th_scale*normcdf(type1_3(3,tt)+sex*sig_th1+sim_factors*sig_th2_n)+sig_th_lower_limit;
        theta_i(:,tt)=min(maxima(1), theta_i(:,tt));
        theta_i(:,tt)=max(minima(1), theta_i(:,tt));
        kappa_i(:,tt)=min(maxima(2), kappa_i(:,tt));
        kappa_i(:,tt)=max(minima(2), kappa_i(:,tt));
        sig_th_i=min(maxima(3), sig_th_i);
        sig_th_i=max(minima(3), sig_th_i);    
        
        
        %%% compares an individual's level of risk aversion given each set
        %%% of simulated factor draws to the threshold level of risk
        %%% aversion proper to a given lottery choice task
        theta_xy_comp=kron(ones(ind_n*draw_n,1), theta_xy);
        theta_diff=(bsxfun(@plus,theta_xy_comp,-theta_i(:,tt)));
        temp=zeros(size(theta_diff));
        for i=1:lottery_n
            temp(:,i)=theta_diff(:,i)./sig_th_i;
        end
        %%% Probability of choosing the risky option i.e. probability that the
        %%% individual's risk aversion coeff is below the threshold for a
        %%% particular lottery. The difference in the coeffs is above
        %%% divided by the individual's sd of the error term so that the
        %%% normal cdf can be used here.
        %%% assumes normal errors on the risk aversion parameter:
        %%% theta~N(0,(sig_th)^2)
        p_risky=normcdf(temp);
        p_safe=1-normcdf(temp);

        % kappa percent of the time and individual makes the wrong choice.
        choice_risky=zeros(size(p_risky));
        for i=1:lottery_n
            choice_risky(:,i)=p_risky(:,i).*(1-kappa_i(:,tt))+p_safe(:,i).*kappa_i(:,tt);
        end
        choice_safe=1-choice_risky;
        sim_l_contrib_lottery=(choice_risky.^(lottery_choice_ind)).*(choice_safe.^(1-lottery_choice_ind));
        
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Time Part %%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% All Types  %%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%
            
    
    
        %%% individual delay aversion based on sex and particular factor draws -
        %%% deterministic part - note that as above, the factor is determined
        %%% by the orthogonal draw AND by X*alpha - a function of observables

        %%% Even though the following parameters are individual-specific - thus
        %%% _i - they vary here within individuals due to the different factor
        %%% draws used in each row for the factor simulation
        r_i=r_scale*normcdf(type4_5(1,tt)+sex*r1+sim_factors*r2_n);

        %%% individual variability of delay aversion based on sex and particular
        %%% factor draws
        sig_r_i=normcdf(type4_5(2,tt)+sex*sig_r1+sim_factors*sig_r2_n)+sig_th_lower_limit;

        r_i=min(maxima(4), r_i);
        r_i=max(minima(4), r_i);
        sig_r_i=min(maxima(5), sig_r_i);
        sig_r_i=max(minima(5), sig_r_i);   

          
    if analytical==1 
     %%%%%%%%%%%% Here is the new alternative approach using an analytical calculation of r_xy thresholds       

                    %%% analytical calculation of indifference thresholds for
                    %%% intertemporal choice tasks

                    %%% first what utility to consider (CRRA, CRRA2, CARA)
            if CRRA2==1
                %%% No longer need to limit risk aversion at 1 as with standard CRRA to keep
                %%% utility positive for discount rate calculations. At
                %%% 10, indifference thresholds are zero 
                theta_limit=min(theta_i(:,tt), 10);   
                theta_limit=max(theta_limit, -.3);            

                %%% utility of either option in period in which it is
                %%% received
                U_1=((x1.^(1-theta_limit)-1)./(1-theta_limit));  
                U_2=((x2.^(1-theta_limit)-1)./(1-theta_limit));

                %%% log case only relevant for John's U as otherwise restrict
                %%% theta to <1
                %%% for theta=1, replaces the above with a ln(x1)

                %%% as with standrd CRRA, at theta=1, U=log(X1)
                [xxx,yyy]=find(theta_limit==1);
                U_1(xxx,:)=repmat(log(x1),size(xxx,1),1);  
                U_2(xxx,:)=repmat(log(x2),size(xxx,1),1);   
            elseif CARA==1
                %%% as above, relaxes upper limit on theta for discount
                %%% rate calculation as CARA stays 
                %%% positive for positive $ amounts.
                theta_limit=min(theta_i(:,tt), 1);   
                theta_limit=max(theta_limit, -.01);            

                U_1=(1-exp(-theta_limit.*x1))./theta_limit;   
                U_2=(1-exp(-theta_limit.*x2))./theta_limit;

                %%% for CARA, at theta=0, U=x1
                [xxx,yyy]=find(theta_limit==0);
                U_1(xxx,:)=repmat(x1,size(xxx,1),1);  
                U_2(xxx,:)=repmat(x2,size(xxx,1),1);               
            else
                %%% Need to keep utility positive, so theta<1
                theta_limit=min(theta_i(:,tt), .99);
                theta_limit=max(theta_limit, -.3);            
                U_1=((x1.^(1-theta_limit))./(1-theta_limit));
                U_2=((x2.^(1-theta_limit))./(1-theta_limit));
            end

            %%% second, what discounting to consider: exponential/hyperbolic
            if hyper==1
                r_xy_comp=((U_2./U_1)-1)./(t2-t1);    
            else
                r_xy_comp=(U_2./U_1).^(1./(t2-t1))-1;
            end
 

    end



        %%% Assumes lognormal errors on the delay aversion parameter, so log:
        %%% ln(r/sig_r)~N(0,1).   
        r_i2=r_i.^2;
        %%% variance of the error term (with a lognormal distrib) on the discount
        %%% rate
        var=(sig_r_i).^2;
        %%% Formulas for mean and sd of normal distrib of which the lognormal
        %%% distrib with mean=r_i and sd=sig_r is the exponent. Done to get pdf,
        %%% cdf of the lognormal distrib which in Matlab takes as parameters the
        %%% mean and sd of the corresponding normal distrib of which it is the log.
        % mean of standard normal distrib
        mu=log(r_i2./sqrt(var+r_i2));
        % sd of standard normal distrib. Adds the small number at the end
        % for purposes of numerical optimization
        sigma=sqrt(log(var./r_i2+1))+(1*10^-20);

        mu_full_time=repmat(mu,1,time_choice_n);
        sigma_full_time=repmat(sigma,1,time_choice_n);

        %%% Probability of choosing the later option i.e. probability that the
        %%% individual's discount rate is below the threshold for a
        %%% particular choice task, given the individual's coefficient of
        %%% risk aversion. Need to take the log here of the 
        %%% resulting difference divided by sig_r as it is the log of the
        %%% error which is normally distributed. 
        p_distant=logncdf(r_xy_comp,mu_full_time,sigma_full_time);
        p_near=1-p_distant;

        %%% If have an extra intercept for time mistakes, need a
        %%% recalculation of Kappa here
        if K_time==1
            kappa_i(:,tt)=0.5*normcdf(kappa_extra+type1_3(2,tt)+sex*kappa1+sim_factors*kappa2_n);

        end
        
        % kappa percent of the time and individual makes the wrong choice.
        choice_distant=zeros(size(p_distant));
        for i=1:time_choice_n
            choice_distant(:,i)=p_distant(:,i).*(1-kappa_i(:,tt))+p_near(:,i).*kappa_i(:,tt);
        end
        choice_near=1-choice_distant;
        sim_l_contrib_time=(choice_distant.^time_choice_ind).*(choice_near.^(1-time_choice_ind));   

        
        %%%%%%%%%%%%%%%%%%%%%% Combining the likelihood contributions from the various
        %%%%%%%%%%%%%%%%%%%%%% types of measures and lottery and
        %%%%%%%%%%%%%%%%%%%%%% intertemporal task choices 


       sim_l_contrib_total=[sim_l_contrib_all_factors sim_l_contrib_lottery sim_l_contrib_time];
       % for numerical estimation purposes

       sim_l_contrib_total_plus=sim_l_contrib_total+1*10^(-20);

        % Individual likelihood contribution conditional on factor draws
        sim_ind_draw_contrib=prod(sim_l_contrib_total_plus,2);
        
        % individual likelihood of measures and choices averaged
        % over simulated factor draws
        av_l_contrib=cell2mat(accumarray(ind_id,sim_ind_draw_contrib,[],@(x){sum(x)}))/draw_n;

        % In order to avoid zeros for the log
        av_l_part(:,tt)=av_l_contrib+1*10^(-200); 
        sim_l_contrib_lottery=NaN(ind_n,lottery_n);
        sim_l_contrib_time=NaN(ind_n,time_choice_n);     
        clear sim_l_contrib_total sim_l_contrib_total_plus sim_ind_draw_contrib av_l_contrib sig_th_i theta_xy_comp theta_diff p_risky p_safe choice_risky choice_safe r_i sig_r_i theta_i_cat r_xy_comp r_i2 var mu sigma mu_full_time sigma_full_time p_distant p_near choice_distant choice_near

    end

    %%%%%%%%%%%%%%%%%%%%%%% Combining the likelihood contributions of the
    %%%%%%%%%%%%%%%%%%%%%%% various unobserved types and weigthing them by each type's
    %%%%%%%%%%%%%%%%%%%%%%% estimated population prevalence
    
    av_l=zeros(ind_n,1);
    for tt=1:n_type
        av_l=av_l+p_type(tt)*av_l_part(:,tt);
    end  
    
    log_likelihood=-nanmean(log(av_l));
    
    %%% will capture the % of nan and inf in the final average likelihoods
    nan_pct=1-mean(isfinite(av_l));
    if nan_pct>.05
        log_likelihood=NaN;
    end
    
    %%%%%%%%%%%%%%%%%%%%%% Below code gets Parameter Values for the Average
    %%%%%%%%%%%%%%%%%%%%%% Person (Table 3 of paper)
    
    mean_sim_factors=mean(sim_factors);  
    mean_f_sim_factors=mean(sim_factors(sex==1));
    mean_m_sim_factors=mean(sim_factors(sex==0));
    
    %%% With types, mean intercept values are taken as a weighted average
    %%% across types
    
    theta_mean=0;
    kappa_mean=0;
    sig_th_mean=0;
    r_mean=0;
    sig_r_mean=0;
    
    if K_time==1
        av_sim_type_prefs=NaN(n_type,6);
    else
        av_sim_type_prefs=NaN(n_type,5);
    end
        
    for tt=1:n_type
        theta_mean=theta_mean+type1_3(1,tt)*p_type(tt);
        kappa_mean=kappa_mean+type1_3(2,tt)*p_type(tt);
        sig_th_mean=sig_th_mean+type1_3(3,tt)*p_type(tt);
        r_mean=r_mean+type4_5(1,tt)*p_type(tt);
        sig_r_mean=sig_r_mean+type4_5(2,tt)*p_type(tt);
        
        %%% for robustness with a separate mistake parameter for temporal
        %%% tasks, gets two mistake parameters (risk and time
        if K_time==1
            av_sim_type_prefs(tt,:)=[ [(type1_3(1,tt)+theta1*mean(sex)+sum(theta2_n.*mean_sim_factors'))] sig_th_scale*normcdf([(type1_3(3,tt)+sig_th1*mean(sex)+sum(sig_th2_n.*mean_sim_factors'))]) r_scale*normcdf([(type4_5(1,tt)+r1*mean(sex)+sum(r2_n.*mean_sim_factors'))]) normcdf([(type4_5(2,tt)+sig_r1*mean(sex)+sum(sig_r2_n.*mean_sim_factors'))]) 0.5*normcdf([(type1_3(2,tt)+kappa1*mean(sex)+sum(kappa2_n.*mean_sim_factors'))]) 0.5*normcdf([(kappa_extra + type1_3(2,tt)+kappa1*mean(sex)+sum(kappa2_n.*mean_sim_factors'))]) ];
        else
            av_sim_type_prefs(tt,:)=[ [(type1_3(1,tt)+theta1*mean(sex)+sum(theta2_n.*mean_sim_factors'))] sig_th_scale*normcdf([(type1_3(3,tt)+sig_th1*mean(sex)+sum(sig_th2_n.*mean_sim_factors'))]) r_scale*normcdf([(type4_5(1,tt)+r1*mean(sex)+sum(r2_n.*mean_sim_factors'))]) normcdf([(type4_5(2,tt)+sig_r1*mean(sex)+sum(sig_r2_n.*mean_sim_factors'))]) 0.5*normcdf([(type1_3(2,tt)+kappa1*mean(sex)+sum(kappa2_n.*mean_sim_factors'))]) ]; 
        end
        
    end
        
    if K_time==1
        av_sim_person_prefs=[ [(theta_mean+theta1*mean(sex)+sum(theta2_n.*mean_sim_factors'))] sig_th_scale*normcdf([(sig_th_mean+sig_th1*mean(sex)+sum(sig_th2_n.*mean_sim_factors'))]) r_scale*normcdf([(r_mean+r1*mean(sex)+sum(r2_n.*mean_sim_factors'))]) normcdf([(sig_r_mean+sig_r1*mean(sex)+sum(sig_r2_n.*mean_sim_factors'))]) 0.5*normcdf([(kappa_mean+kappa1*mean(sex)+sum(kappa2_n.*mean_sim_factors'))]) 0.5*normcdf([(kappa_extra+kappa_mean+kappa1*mean(sex)+sum(kappa2_n.*mean_sim_factors'))])  ]; 
        av_sim_women_prefs=[ [(theta_mean+theta1*1+sum(theta2_n.*mean_f_sim_factors'))] sig_th_scale*normcdf([(sig_th_mean+sig_th1*1+sum(sig_th2_n.*mean_f_sim_factors'))]) r_scale*normcdf([(r_mean+r1*1+sum(r2_n.*mean_f_sim_factors'))]) normcdf([(sig_r_mean+sig_r1*1+sum(sig_r2_n.*mean_f_sim_factors'))]) 0.5*normcdf([(kappa_mean+kappa1*1+sum(kappa2_n.*mean_f_sim_factors'))]) 0.5*normcdf([(kappa_extra+kappa_mean+kappa1*1+sum(kappa2_n.*mean_sim_factors'))])  ]; 
        av_sim_men_prefs=[ [(theta_mean+theta1*0+sum(theta2_n.*mean_m_sim_factors'))] sig_th_scale*normcdf([(sig_th_mean+sig_th1*0+sum(sig_th2_n.*mean_m_sim_factors'))]) r_scale*normcdf([(r_mean+r1*0+sum(r2_n.*mean_m_sim_factors'))]) normcdf([(sig_r_mean+sig_r1*0+sum(sig_r2_n.*mean_m_sim_factors'))]) 0.5*normcdf([(kappa_mean+kappa1*0+sum(kappa2_n.*mean_m_sim_factors'))]) 0.5*normcdf([(kappa_extra+kappa_mean+kappa1*0+sum(kappa2_n.*mean_sim_factors'))])  ]; 
    else
        av_sim_person_prefs=[ [(theta_mean+theta1*mean(sex)+sum(theta2_n.*mean_sim_factors'))] sig_th_scale*normcdf([(sig_th_mean+sig_th1*mean(sex)+sum(sig_th2_n.*mean_sim_factors'))]) r_scale*normcdf([(r_mean+r1*mean(sex)+sum(r2_n.*mean_sim_factors'))]) normcdf([(sig_r_mean+sig_r1*mean(sex)+sum(sig_r2_n.*mean_sim_factors'))]) 0.5*normcdf([(kappa_mean+kappa1*mean(sex)+sum(kappa2_n.*mean_sim_factors'))]) ]; 
        av_sim_women_prefs=[ [(theta_mean+theta1*1+sum(theta2_n.*mean_f_sim_factors'))] sig_th_scale*normcdf([(sig_th_mean+sig_th1*1+sum(sig_th2_n.*mean_f_sim_factors'))]) r_scale*normcdf([(r_mean+r1*1+sum(r2_n.*mean_f_sim_factors'))]) normcdf([(sig_r_mean+sig_r1*1+sum(sig_r2_n.*mean_f_sim_factors'))]) 0.5*normcdf([(kappa_mean+kappa1*1+sum(kappa2_n.*mean_f_sim_factors'))]) ]; 
        av_sim_men_prefs=[ [(theta_mean+theta1*0+sum(theta2_n.*mean_m_sim_factors'))] sig_th_scale*normcdf([(sig_th_mean+sig_th1*0+sum(sig_th2_n.*mean_m_sim_factors'))]) r_scale*normcdf([(r_mean+r1*0+sum(r2_n.*mean_m_sim_factors'))]) normcdf([(sig_r_mean+sig_r1*0+sum(sig_r2_n.*mean_m_sim_factors'))]) 0.5*normcdf([(kappa_mean+kappa1*0+sum(kappa2_n.*mean_m_sim_factors'))]) ]; 
    end
    
    av_sim_prefs=[av_sim_person_prefs;av_sim_women_prefs;av_sim_men_prefs;av_sim_type_prefs];
    
  
    

    
end